{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Finite Element Methods\n",
    "\n",
    "We use the linear finite element method for solving the Poisson equation as an example to explain the main ingredients of finite element methods. We recommend to read \n",
    "\n",
    "- [Introduction to Finite Element Methods](http://www.math.uci.edu/~chenlong/226/Ch2FEM.pdf)\n",
    "- [Progamming of Finite Element Methods](http://www.math.uci.edu/~chenlong/226/Ch3FEMCode.pdf)\n",
    "\n",
    "and do the following project\n",
    "- [Projects: Linear Finite Element Methods](../project/projectFEM.html)\n",
    "\n",
    "The numerical example using the linear element is in\n",
    "- [Linear Element for Poisson Equation in 2D](Poissonfemrate.html)\n",
    "\n",
    "and more elements and equations can be found\n",
    "- [List of Examples](femexamplelist.html)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Variational formulation\n",
    "\n",
    "The classic formulation of the Poisson equation reads as\n",
    "\n",
    "$$ - \\Delta u = f  \\text{ in }  \\Omega, \\qquad u  = g_D  \\text{ on }\n",
    "\\Gamma _D,  \\qquad  \\nabla u\\cdot n = g_N  \\text{ on } \\Gamma _N, $$\n",
    "\n",
    "where $\\partial \\Omega = \\Gamma _D\\cup \\Gamma _N$ and $\\Gamma _D\\cap \\Gamma _N=\\emptyset$. \n",
    "We assume $\\Gamma _D$ is closed and $\\Gamma _N$ open.\n",
    "\n",
    "Denoted by $H_{g_D}^1(\\Omega)=\\{v\\in L^2(\\Omega), \\nabla v\\in L^2(\\Omega) \\text{ and } v|_{\\Gamma _D} = g_D\\}$. Multiplying the Poisson equation by a test function $v\\in H_{0_D}^1$ and using integration by parts, we obtain the weak formulation of the Poisson equation: find $u\\in H_{g_D}^1(\\Omega)$ such that for all $v\\in H_{0_D}^1$:\n",
    "\n",
    "$$ \n",
    "a(u,v) := \\int _{\\Omega} \\nabla u\\cdot \\nabla v\\, {\\rm dxdy} = \\int _{\\Omega} fv \\, {\\rm dxdy} + \\int _{\\Gamma _N} g_N v \\,{\\rm d}S.\n",
    "$$\n",
    "\n",
    "Let $\\mathcal T$ be a triangulation of $\\Omega$. We define the linear\n",
    "finite element space on $\\mathcal T$ as \n",
    "\n",
    "$$\n",
    "\\mathcal V_{\\mathcal T} = \\{v\\in C(\\bar \\Omega) : v|_{\\tau}\\in \\mathcal P_k, \\forall \\tau \\in \\mathcal T\\}. \n",
    "$$\n",
    "\n",
    "where $\\mathcal P_k$ is the polynomial space with degree $\\leq k$. \n",
    "\n",
    "The finite element method for solving the Poisson\n",
    "equation is to find $u\\in \\mathcal V_{\\mathcal T}\\cap H_{g_D}^1(\\Omega)$ \n",
    "such that for all $v\\in \\mathcal V_{\\mathcal T}\\cap H_{0_D}^1(\\Omega)$:\n",
    "\n",
    "$$\n",
    "a(u,v) = \\int _{\\Omega} fv \\, {\\rm dxdy} + \\int _{\\Gamma _N} g_N v \\,{\\rm d}S.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Finite element space\n",
    "\n",
    "We take linear finite element spaces as an example. For each vertex $v_i$ of $\\mathcal T$, let $\\phi _i$ be the piecewise linear function such that $\\phi _i(v_i)=1$ and $\\phi _i(v_j)=0$ when $j\\neq i$. The basis function in 1-D and 2-D is illustrated below. It is also called hat function named after the shape of its graph.\n",
    "\n",
    "Then it is easy to see $\\mathcal V_{\\mathcal T}$ is spanned by $\\{\\phi\n",
    "_i\\}_{i=1}^{N}$ and thus for a finite element function $v=\\sum\n",
    "_{i=1}^Nv_i\\phi _i$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA\nB3RJTUUH4gkHFSof/xHqpwAAACR0RVh0U29mdHdhcmUATUFUTEFCLCBUaGUgTWF0aFdvcmtzLCBJ\nbmMuPFjdGAAAACJ0RVh0Q3JlYXRpb24gVGltZQAwNy1TZXAtMjAxOCAxNDo0MjozMROMmUQAACAA\nSURBVHic7d1/dFNlnj/wJzS0oSNI4lTa4kpqBzOox36lxkq1mkQt6oHRmXMc45kdKaOzYsVhUUR3\nYLapq8c94LoUBFnO2Wm7c6bW+boOQ5cRe4S2OEi/LcU9M1Z+1JiLdFooNQlF2nRayPePp1xCkiZp\nknvvc5/7fh3/oO01eSjPzec+Pz6fRxcMBgkAAIDSpindAAAAAEIQkAAAgBEISAAAwAQEJAAAYAIC\nEgAAMAEBCQAAmICABAAATEBAAgAAJiAgAQAAExCQAACACQhIAADABAQkAABgAgISAAAwAQEJAACY\ngIAEAABMQEACAAAmICABAAATEJAAAIAJCEgAAMAEBCQAAGACAhIAADABAQkAAJiAgAQAAExAQAIA\nACYgIAEAABMQkAAAgAnqDkjffvvtyZMnlW4FgFTQw0FT1B2Q3n777a1btyrdCgCpoIeDpqg1INXU\n1DzxxBO1tbVKNwRAEujhoEF6pRuQpNLS0uLi4t27dweDQaXbApB+6OGgQWoNSFarlRDy+eefC4Kg\ndFsIIUTwBgghZpNB6YYAJ1jr4QAyUGtASsRPf/rTjo4OQsjKlSuff/55Sd+rvrPf1expqbzNVmiU\n9I0ARDL0cMEbqO/srzvUH/rNitvzqhYXSPF2oHE8B6SOjo5jx47J816CL0DoOKlQnjcEkLyHt7p9\nyxuPCN5AhTUv9PuuZk/dof5a5wI8fkF68RyQ5HfCG1C6CQDpUf2Rx9XsMZsMkeP+qvIC+zuHlzce\nwVAJ0gsBKT1a3T5yaZwEoHY0GtkKZ7dULoz8qdlkaHl2IZ2mFnyBWucC+VsIXFLrtm82Cd4RpZsA\nkAYxohFlNhmqFhfYCmfTRzGAtFD3CGnFihVKN2EC3WWHERKklyI9vK6znxCy7Mp1o6iqFhfYt31W\n19lfkcDFAHFhhJQGApaOgCP1nf2EkERijNk4Q7weIHUISGkg+CZm6hCZgAOtbn+CIx6zyWArnI2J\nAUgXBKQ0oHHIVjibICaBygnegNlkuLdwdtj3W90+3Yv7ljceCfs+ndlDt4e0QEBKG7o1VhwtAagU\nTYYN+2bbl35CSF1nf9guhja3n8Yw+doH/EJASoM2t58QMs9kIHhUBJWbbBZO/E7YIAk7GiCNEJDS\nho6QkBsLarfMmid4A2EjoVa3r8KaV+tcIHgD1R956DfpfrzI+T2A5CAgpYHgHcGUBXCDPlqJUYcQ\n0ur2Cd7AvYWzbYVGW+HsukP9dCaAzg2ggBCkCwJSGgi+gK3QaDYZzCYD8gRB7cwmg6u8oNXtF2fn\nLm3bMdJ8WMEbqG722Lcdruvsd5UX4GkM0kXdibGMEGuqmo24M4EHdO+cq9nT6vbVOhfQkRANPLZC\nY4U1j07WucoLUMsO0ggBKT3EUIScDOAAHQkts+bRIqp0H504YKLRiF6jaDOBN5iySxWdo6Nb7Mym\nGdhlB9ygRVSryieiTqvbR/9zlReE7W4ASAuMkFIVGoHoOAlpGcANs8lA3IQQUlVeELq9W/AG6O6G\nZdY89HZIF4yQ0oNuNJpIRUJuLHAk6lY6s8lQ67yJ7m5QqF3AIQSkVIUmHuFREfhDJ6Uj+7bZZKC7\nG7CzFNIFASlVdBcDvV1p8WNaZAWAD5FHmIuqygtCNzsApAgBKVXIigWOxa7FYDYZqsoLsLsB0gUB\nKVWCLyDu+aaRCTMYwI0Tl1JiJ7ugwpoXWrsBIBUISKkSvAGzaYb4pQ11vYAjoTPSk6G7G5Y3fiFX\no4BbCEhpEFagAbmxwA1aUzX2NWKpIcwNQIoQkFISmhVLITcWuCHWVI17Jc1Gwu4GSBECUkoiY4+Y\nG6tEcwDSSYi3gCQSdzcgJkEqEJDSIPSORW4scCO0pmpcdHcDrS0kcbuAWwhIKYk8jg9bwIEbiSwg\nhZqo3YAt4JAsdgNSIBDo6ekZGhqa7ILh4eGenp7h4WE5WxUmcg8ScmMhQYn08OPHj8e4QGpTnXkW\ndzfQ7CWAqWI0IDU1NZWVla1Zs8Zut2/cuDHygm3bttlstrVr195zzz2//e1v5W8hhaxYSE7cHv7+\n++/bbLaXX365rKzs9ddfl7+FyR1PTnc3oMAdJCnIHp/Pd+utt7a3tweDwYGBgTvuuOPAgQOhF/zp\nT3+68847z5w5EwwG3W73woULv/zyy8jXufHGG6Vuqvm1A7atXWHfJC/sjfwmgChuDz9//vyCBQsO\nHz4cDAb7+/tvvvlm+ucwkvbwine/IC/s9XwzMtX/sbajj7ywt+LdL6RoFfCNxRFSR0dHbm5uSUkJ\nISQnJ8fhcOzfvz/0gq6urkWLFn33u98lhNxwww233nrrJ598okhTw7JiKYyZILa4PVyn0+n1+jlz\n5hBCZs2aNX369MzMTJkbOVlN1biwuwGSxmJAOn36dF7e5aXU3NzcgYGB0Atmz57d19dH/xwMBk+d\nOiV+Kb/IY8vNRgNyYyGGuD18xowZa9eurays3LZt21NPPbV06dKbbrpJ5kbGqKkaF3Y3QHJYDEjj\n4+MZGRnil3q9fmxsLPSCBx988OjRo2+99dZnn3326quvnj9/fnR0NOpLWSwWi8WyZcsWKdpJl3zn\nRSnLHz5mAggVt4dfuHDh66+/Pnv2bF9f3+joaG9vr9frjfpSEvXw5BaQRNjdAMlh8cTYrKysQODy\nCGNkZCQrKyv0gmuvvfbdd9/9j//4j7feeuuBBx5wOBwmkynqSx07dky6dk42I2E2GgRvAOfGwmTi\n9vCurq7du3fv2bNn5syZwWBw2bJl7733XmVlZeRLSdTD49ZUjWuZNa/uUH91syfpYRZoEIsjpLlz\n5wqCIH7p8Xiuu+660Av8fv/Fixc3bdr0m9/85sknnzx69GhhYaHcrbwk8qZFbizEFreH9/T05Ofn\nz5w5kxCi0+mKi4vdbrecLUykpmpsZpOh1rkAtRtgSlgMSFardXR0tLGxkRDS3d3d1tZmt9sJIQ0N\nDV1dXYSQjIyMn/3sZ0eOHCGE/OlPf+rp6XE4HPK3MzIrlsLACGKL28OLior+8pe/HDhwgBDi9Xo/\n/PDD4uJiOVtY19mf+sjGVmjE7gaYEhan7LKzszds2PDLX/5y8+bN586dW7lyZVFRESFk8+bNjz/+\neHFx8cyZM1966aVVq1aNjo5evHhx27ZtM2YosGwz2VMkzY0VvAGi2LANmBa3h99yyy3r1q1btWqV\nwWAYHh5+5JFHnE6nbM2j8SPpBaRQtc6bCl7/tPojj60y+dk/0A5dMBhUug3RBYPBwcFBo9Go108a\nNX0+n9E4aUe3WCySriHZtx0WfAHPutKw7wveQMHrn7rKC6oWF0j37qB2cXv4xYsX6QXTp0+PeoFE\nPbyus3954xHPutK0jPXpq9U6F2AxCeJiccqO0ul0OTk5MaIRISRGNJJB6FmxocwTa0jY+Q2xxO3h\n06ZNu/baayeLRtKZUk3VuOjEHWo3QCLYDUiqMNkOb7PJIHixqQFUaao1VWMzmwxVi3EyBSQEASl5\nMUpPIjcW1Cvtp3lhdwMkCAEpSfSmnWzhF7mxoFIppsROhtZuwCAJYkNASlLsZz0xN1a29gCkBV1A\nSiUlNioxLQm1GyAGBKSUxF74RW4sqE7SNVXjEnc34EENJoOAlCSaFUtTjiLd+700z3gAyCOVmqqx\nibsbsOMOJoOAlKREaqvgSRDURaIFJJGt0Fhhzavr7MfuBogKAUkSdOQ0WW0hADalXlM1rqryAkII\ndjdAVAhISWp1+2yTP0giNxbUKPWaqnGJuxtwWhJEQkBKXuy93ciNBdVJS03VuOjuhrpD/ZjThjAI\nSEmKey8hNxbUJY01VWPD7gaYDAJSMoSJLXY4ZgL4IUi/gCTC7gaICgEpGTTBKPLw8lC2QiNmJEBF\n0ltTNS7sboBICEjJmBghJXDr4gEQ1CK9NVXjwu4GiISAlIzYWbEUcmNBXeQf0GN3A4RBQEpG4rtj\ncaeBKkidEhuV2WSgRVexuwEoBCSpIDcWVESimqpxmU0G7G4AEQJSMmJnxVLIjQUVaXX7zCaDbDsa\nQlWVF5hNBuxuAIKAlLRETjxCbiyoheANyD88oswmQ1V5AXY3AEFASk6CK0PIjQVVUGQBKVSFNQ+7\nG4AgICUBWbHAGRlqqsZ16UjZLxRsAyhOr3QD1CeRrFjKVmh0YfsQME+GmqpxmU0GV3mBq9lT19l/\nwhuoO3TFwbIVt+dVLS5Qqm0gG3YDUiAQOHny5Jw5c2bNmhX1guHh4d7e3tzc3MkukEjiWbHi9cre\n6sCmuD2cENLf33/11VdnZ2dL2hJ5aqrGtcyaV3eon+5uCGuPq9lTd6i/1rlA2WEcSI3RgNTU1PTq\nq6/m5+f39vY6nc6XXnop7IL3339/w4YNc+fO/eqrr3784x+vW7dO5hbGzoql7v3ebNJMBN8IAhKE\nidvDT5w4sWLFigsXLpw/f95ut7/22msStUS2mqpx1XdOrCFVWPNqnQtCf1RVXmB/5/DyxiMYKnEu\nyB6fz3frrbe2t7cHg8GBgYE77rjjwIEDoRecP39+wYIFhw8fDgaD/f39N998M/1zmBtvvFGK5lW8\n+wV5YW8iV7Z86SUv7K3t6JOiGaBecXv4xYsXH3jggV27dgWDwZGRkfvuu6+zszPyddLSw2s7+sgL\nez3fjKT+Uqlw7fmKvLDXtrXLtrXL/NqBli+9YRd4vhmh11S8+4UiLQQZsLipoaOjIzc3t6SkhBCS\nk5PjcDj2798feoFOp9Pr9XPmzCGEzJo1a/r06ZmZmcq0NSbkxkJUcXv44cOHdTrd0qVLg8GgwWBo\nbm4uLi6WqDEy11SdjKvZYyuc3VK5cKJ2Q8QWcHpoha1wNlJoOcZiQDp9+nRe3uUZ5Nzc3IGBgdAL\nZsyYsXbt2srKym3btj311FNLly696aabZGseTSFM/Hrs/IYwcXv48ePHLRaLy+UqKSm58847t2zZ\nIl1jZK6pGhXdd77Mmkcu7W5odfvpN8PQg5Si/gg4wGJAGh8fz8jIEL/U6/VjY2OhF1y4cOHrr78+\ne/ZsX1/f6Ohob2+v1+uN+lIWi8VisaT9fk5wZZWmviM3FsLE7eF+v3/fvn25ubmffPJJbW3t7373\nu127dkV9qdR7OAupP/Wd/SRkI8Mya57ZZIha4I7OOtQjIHGKxU0NWVlZgcDlm2RkZCQrKyv0gq6u\nrt27d+/Zs2fmzJnBYHDZsmXvvfdeZWVl5EsdO3Ys7c0TvAFSmOjFSFeCSHF7uMFguPrqq5955hmd\nTrdgwYLHHnuspaXlkUceiXypFHu44imxVKvbHzpKo7UbljceWd54JGx3g9lksBXOxqwDr1gcIc2d\nO1cQBPFLj8dz3XXXhV7Q09OTn58/c+ZMQohOpysuLna73fK0LYmsWNw8ECZuDzebzZmZmTqdjn6Z\nmZl54cIFKVqiVE3VUDQvIiwo0toNUZeL6MweCwM7SDsWA5LVah0dHW1sbCSEdHd3t7W12e12QkhD\nQ0NXVxchpKio6C9/+cuBAwcIIV6v98MPP5RuyTdM4lmxlNk0A3cOhInbwxctWjQ8PPzhhx8SQrxe\n786dO202mxQtoWWCFd/RIHgDkbNw9MzlyJjU5vYjt49XLAak7OzsDRs2bNq0qbS01Ol0rly5sqio\niBCyefNmuhnplltuWbdu3apVq+6+++77779/0aJFTqdTnrZNNSuWjqUQkyBU3B5uMBi2b99eU1Nj\nt9sfeughu93+6KOPStESwRtIpEywpCabhaNHXEbeO4yk8YIUWFxDIoQ4HI6DBw8ODg4ajUa9fqKR\n7e3t4gU//elPf/KTn9ALpk+fLnPzEsmKpeZNHEKB3Fi4Qtweftttt+3Zs8fv98+cOTN0B0QaMbKA\nRAhZZs1b3nik1e0LnTyMmjXBTptBCiyOkCidTpeTkyPeq5GmTZt27bXXyhyN6Jx74iZORcIICSLE\n7eGEkNmzZ0sUjQgbNVUp2oaw3CM6cgqbsmNh0Qukw25AYtlUpuyQGwuMop/1LIzdxdyjyGP6xKk8\nwRuwbztc19nvKi9goc0gBUan7JgleDH5BpwI22ytLLp3ztXsaXX7xCKqtF4+DZzLG48I3oCrvAC1\n7DiGgDQ1gm9qB2vS3NhWt6+K4C4ChrBTU5WilYGWWfNoEVV6l9Gk8vrOU3Wd/WaToaXyNkzW8Q0B\naWqmlBVLITcWGCQws4AUymwytDy7sNXtCy3T0Or2YWCkEQhIU5ZEgEFuLLCGkZqqkcwmQ4UpT5xL\n1L24z2w0IBppBDY1TAGd5Ug8K5ZCbiwwiIWaqolAoSBNQUCaguTiCnJjgUGCN6CKyWRarwG3j0Yg\nIE3ZVKfdxdxYaZoDMGU0vXSqY31FTNRrwO2jDQhIU5BcOhGD0/SgcSpKL6WZfG1fTi0hHVQKAWkK\n6Fz2VAMM7ihgDSM1VROHU2I1AgFpCpAVC3xgoaZqgnAAkqYgIE2B4EtmHZjGMDziASNUV59ULbET\nUoeANAVJP1fa1HPzA/fYqamaILPREPVgJOAPAtLUJL1TFnMOwAh2aqomaLKDkYA/CEiJSi4rlkJu\nLLCDqZqqiUDJfO1AQEpUKhEFubHACNZqqiYi6sFIwCUEpKlJbuYdubHACDZrqiYCk95agICUqFRm\nDFQ0Xw98Y7amamy0gBAGSdxDQEpUclmxFHJjgRFqqakahu5rAO4hICUKWbHAAbXUVA2DRzqNQEBK\nVHJZsRRyY4EFKqqpGgZ3kEYgICUqxWorGF2B4ugCkhqn7AgORtIGBCSZmI0G3E6gLFpTVelWJInu\na1C6FSAtdgNSIBDo6ekZGhqK/NH4+Pi5KwUC0vZUeiekkr2BelwQJkYPD3Xq1KnBwcG0vKOKaqpG\novsaMGvHN0YDUlNTU1lZ2Zo1a+x2+8aNG8N++vHHH5eFKCkp+Zd/+RdJ25P6bUDrceERD6jYPVx0\n7ty5J554oqmpKfV3VF1N1TDY16AFeqUbEIXf71+/fv2OHTtKSkrOnDmzZMmSu+66q7S0VLzgwQcf\nfPDBB+mfBwcHn3766aefflqGhqWSTijmxmIxCeL2cFF1dfVVV12VljdV0aF8MbS6fVWkQOlWgFRY\nHCF1dHTk5uaWlJQQQnJychwOx/79+6NeeeHChRUrVrzwwgsFBdL20dTraCEOgSjBHt7U1GQ0GouL\ni9PypoJ3hKi5H+JgJC1gMSCdPn06L+/yRqDc3NyBgYGoV/7ud7+bOXPmPffcM9lLWSwWi8WyZcuW\nFJuUSlbsFa+DKTtIrIf39fXV19e/+OKLsV8q8R6uupqqkdS7AAYJYnHKbnx8PCMjQ/xSr9ePjY1F\nXhYIBDZv3rx169YYL3Xs2LG0NCn1rFhULAZR3B5+8eLFV155Zd26dQZDnF6XYA9X+wISJR6MpPaJ\nR5gMiyOkrKys0F1zIyMjWVlZkZft27cvOzt74cKFMjQplaxYyjyxhoSABPF7eH19vdFo/M53vnP8\n+HGfz3fmzJm+vr7U31ftn+M4GIl7LI6Q5s6dKwiC+KXH45k/f37kZU1NTeLWBhmkPl1gNhnoPD5o\nXNwe7vf7PR7PmjVrCCGnTp3KzMwcGxtbt25d0u+o0pqqYTDNwD0WR0hWq3V0dLSxsZEQ0t3d3dbW\nZrfbCSENDQ1dXV3iZX/+85+LiorkaVJaHsqQGwtU3B6+evXqXZc8/PDDTz31VCrRiKi2pmoYs8lg\nNhmQisQxFgNSdnb2hg0bNm3aVFpa6nQ6V65cSQPP5s2bxc1IZ8+eHRwcLCwslKE9qWfFUliSBSqR\nHp5eKq2pGglPdXxjccqOEOJwOA4ePDg4OGg0GvX6iUa2t7eLF1x99dXp2rAQV7qeyMTcWLXPnEDq\n4vZwkcvlSvG91FtTNZKt0Ohq9uAm4hWLIyRKp9Pl5OSI96ri0nUD4NxYoGTr4aquqRpmYl8DbiJO\nsRuQ2EEXUemCaipwyBgoQtU1VcOggBDfEJDiS1dW7MSrYY8QyEvVNVXD4GAkviEgyQebVkF+fKTE\nhkIBIY4hIMWXrhkP5MaC/PioqRoKByNxDAEpIema8UBuLMhM7TVVI+FgJI4hIMWXxscxZFGAzDio\nqRoG+xo4hoAUhzCxxY6fB0zQDv4WkEQYIXEJASkOmvGQrqRCTH+D/HhaQCI4GIlrCEhxTIyQ0joF\nj4c7kAcfNVVBOxCQ4khXViyF3FiQEx81VSPRmQY82PEHASmO9GbFTrwmZu1AFtzUVA2DBzteISDJ\nCrmxIBueaqqGwUY7XiEgxZHeOmDIjQXZ8FRTNQwORuIVAlJ86a0DhtxYkAdPNVUjIaWPSwhIcaR9\nvQc3EsiDp5qqkei+BizHcgYBKRZkxYJKcZwSS+FgJC4hIMWS3qxYCrmxIAP+aqqGwb4GLiEgxSJF\nVmzoKwNIhC74c5wSi4ORuISAFF+6smIpTDWADARvgMv9daFQQIg/CEixSFd5BSMkkA73C0gUZr/5\ng4AkN+TGgjw4XkCicDASfxCQYml1+ySahcdUA0hHUzVVsa+BJ+wGpEAg0NPTMzQ0FOOa/v7+4eFh\nSZuR9sdMmmSO3FiI28MTuQWi4rWmahg62YAREk8YDUhNTU1lZWVr1qyx2+0bN26MvODEiRMPPfTQ\nsmXLHnjggfXr10vUDIlmqJHYBHF7eENDQ1lZ2csvv+xwOGpqaqb04rzWVA2Dg5H4o1e6AVH4/f71\n69fv2LGjpKTkzJkzS5Ysueuuu0pLS8ULgsHgz3/+8+eff37p0qWBQGDJkiWHDh26/fbb09sMSbNi\ncRdpWdwe3tPTs3Hjxp07d86bN6+vr+8HP/hBWVnZwoULE3lxjmuqAvdYHCF1dHTk5uaWlJQQQnJy\nchwOx/79+0MvOHz4sE6nW7p0aTAYNBgMzc3NxcXFaW+GFFmxlNk0A7uDtCxuDz969Ghpaem8efMI\nIfn5+ddff/2JEycSfHGOa6pGwsFInGExIJ0+fTov7/LtlJubOzAwEHrB8ePHLRaLy+UqKSm58847\nt2zZMtlLWSwWi8US44IYpMuKpaMuxCTNitvDly5dunXrVvpnQRDcbndRUVHUl4rs4XzXVA2Dg5E4\nw+KU3fj4eEZGhvilXq8fGxsLvcDv9+/bt2/lypWffPLJV1999fTTT5vN5kceeSTypY4dO5ZiY9Kb\nFUvNmziEYkQj+6AgTNweLurq6lq9evWKFStuuOGGqBdE9nDBG+B+w7dILCCknb8y31gcIWVlZQUC\nl0cPIyMjWVlZoRcYDIarr776mWeeycrKWrBgwWOPPdbS0pL2ZtCpDylMnIqEEZJWxe3hhJCxsbE3\n3njjH//xH9evX//ss88m+MoaSYkV4WAkzrA4Qpo7d64gCOKXHo9n/vz5oReYzebMzEydTke/zMzM\nvHDhgkSNkWbKDrmxmha3hxNCnn/+eb1ev3v37lmzZiX+ytzXVI2E81x4wuIIyWq1jo6ONjY2EkK6\nu7vb2trsdjshpKGhoaurixCyaNGi4eHhDz/8kBDi9Xp37txps9nS3gzBiyk1kETcHv7xxx/39vbW\n1NRMKRoRDdRUjYSDkXjCYkDKzs7esGHDpk2bSktLnU7nypUr6Yru5s2b6WYkg8Gwffv2mpoau93+\n0EMP2e32Rx99NO3NEHxSzcVjnkHj4vbwAwcO0I0Mt1zS3NycyCtroaZqGFQr5gmLU3aEEIfDcfDg\nwcHBQaPRqNdPNLK9vV284LbbbtuzZ4/f7585c2bo+nAaCd4AKZTihQlBbqzmxe7hVVVVVVVVU31N\nrS0gUdjXwBMWR0iUTqfLyckR79WoZs+eLVE0oiQNG5j41rhEengStPa5jIOReMJuQFIW7d/Spbsj\nNxbSTlM1VUOhgBA3EJCikzpaIDcW0k4jNVUj4WAkbiAgxSLd7IeYGyvR64MGaaSmaqR5mLXjBQJS\ndFInCWlwXgUkNbGjQZOldOjdhIOROICAFB2dkpYubIhbgyR6fdAaDabEinAwEjcQkKJDViyoi6Zq\nqoahByMp3QpIAwSk6ASftNPx2KsK6SV4A2ZT+gsBqwg22nEAASk6GW5vPNNBumgzJTYUDkbiAwLS\npGTYsIRnOkgLLS8gUdrczcEfBKQopM6KpZAbC+miwZqqYbBLiA8ISFHIEyeQGwvposGaqmFQsJgP\nCEiTknoCBLmxkBZYQKJwMBIHEJCikOfoPC1PsEDaaXkBicLBSBxAQIpC6qxYCrPekBaarakaBgcj\ncQABKQpkxYKKaLamahg84XEAASkKqbNiKeTGQuroJJU2a6qGwQ3FAQSkKGRLesc4DFJEP3+RhUPh\nYCS1Q0BSEvYFQYqQEhsKByOpHQJSONqh5dlEq/HiY5A6LCBFwqydeumVbgBz5OzNZqNhYg0Ac3eQ\nGMEbqO/srzvUTwj5q+O1MW+glfiqP/JULS5QumnKu/d7s0kzafvSjyFjqNA+I6q4PY/BPoOAFJ08\nHVrMjUVAgkS0un3LG4+IdRne6/nzWPY1ZqPB1eypO9Rf61yg8Q9iutEO0+ChwvqMiM0+g4AUTp6s\nWApxCBJX/ZHH1ewxmwwtlbfRD5EPPvhg5BpS67yJEGJ/5/DyxiNsPvbKhh6MJHiRijQhss+IqsoL\nGOwz7AakQCBw8uTJOXPmzJo1K/Knf/vb30ZHR8Uvr7rqKp1Ol5b3lScr9op39AZIoWzvBqoRegvQ\nTxZb4eyWyoXiBSPX3Egu9dWWZxfWd/a7mj2CL1DrXKBYoxmAERIVtc+IzCYDg32G0YDU1NT06quv\n5ufn9/b2Op3Ol156KeyCmpqa9957z2CYCBu///3vc3Jy0vLWcmbF0hkGOcdkoBZht4DrVHHkJ8tY\n9jXiPIzZZKhaXNDq9ml8Sd9WaHQ1e1rdPqZmohQRIxpRDPYZFnfZ+f3+9evXv/3223/4wx/27Nnz\n/vvvf/rpp2HXHD16dNOmTX+6JF3RiMiVFUuZJ9aQEJDgCmG3wPb9HkLIsivX8AYioAAAIABJREFU\nAKLWVK1aXCB4A/RH2oSULIr2gWUJ7MBkqs+wGJA6Ojpyc3NLSkoIITk5OQ6HY//+/WHXHD169Pvf\n/35vb28gkP5Pczl3Y5tNBkx5Q5iwW2Dcch8hJOr27rBxAB1z17Px4aIIFBCiaB9IJCWAqT7DYkA6\nffp0Xt7l32Nubu7AwEDoBd98883Zs2crKiqefPJJq9W6ZcuWNL67zIl1yI2FSGG3wNdB403kZNg1\ndKY3bHp5YlVfwz0KByNRrW5/ZDSK+uHGVJ9hMSCNj49nZGSIX+r1+rGxsdALzp49++CDD/7617/e\nt29fY2Pjf/7nf0bO6VEWi8VisSQeseTMiqWQGwuRQm8BwRsw6cfyxgfCrqlaXHDj/zwb2cPpLI2W\nCxbgIY+mNkZ+jtFtdZHXs9NnWNzUkJWVFToRNzIykpWVFXrBDTfc8Oabb9I/33zzzeXl5e3t7aWl\npZEvdezYsSm9tfwPVjQ3VuY3BcaF3QLe8enH9eaoV0b28Da3X+Op1nRfg8Z/CTQZNnSQRHPw7y2P\n8rTNTp9hcYQ0d+5cQRDELz0ez3XXXRd6QUdHR3Nzs/jlxYsXL1y4kMYGyP8PgxkGCBV6C5hNhvyL\nZ87rZyb4/9Zd+TGkQTgYKeosHP2Qifrhxk6fYTEgWa3W0dHRxsZGQkh3d3dbW5vdbieENDQ0dHV1\nEUICgUB1dXVvby8h5Msvv9y3b999992XlreemJc3yjeNhk1BECnsFtB90ewdz0zkqQXHmRPsayCE\nELLMmid4A6F9ZrI6vEz1GRYDUnZ29oYNGzZt2lRaWup0OleuXFlUVEQI2bx5M91ud8899zz++OOP\nPPLIAw888OSTT65du3bhwkn32k+J/FmxE++LWTsIEXYLPLuklBBS/ZEn7v+I4t8EByMRQi71gdA+\nM1kdXqb6DItrSIQQh8Nx8ODBwcFBo9Go1080sr29XbzgF7/4RWVlpc/nS2MGkiKQGwtRhd0C4x95\nXM2e5Y1HJsuoF7yB5Y1ftLr9rvICFhYDlMXOtjGlmE0GV3lBaJ+JrAjDYJ9hNCARQnQ6Xexgo9fr\n0x6NWt0+m7xDV+TGwmRCbwG6D4rWIIgsiCkW0HSVF7BTl0xBZtOMVremp+zIlX2m4vY8cuW8HJt9\nht2ApBT592EjNxbiolVellnz6M5dGpBOFT1Jd/HWdfZHLaCpWbTYisYLCIX2GVezhxDS5va3XYrT\nbPYZBKQrKFLnFGkTkCBaELPV7atu9hBCRq65ka6UMPWQywIcjCSifcb+zuGwPQ5s9hkEpMuEiS12\nTMylAkRlNhkqTHl0ddpisUw1004j5Nwoyz6zyUDPQ2KkpHcMLO6yUwpNXJgn++KerdCIXXYAaUQT\ncTS+0U7E1Mbu2BCQLhOiFQeTDW4egPTCTDhFN/GqYvYSAeky+bNiKeTGAqQdnXjAcx5RLr0yCQhI\nlyn7z4ZZO4A0wnOeaLKUWAYhICkPubEAaYcCQlSr2yd4A6pYQCIISKHkz4qlkBsLkHY4GIkS1LOA\nRBCQwih1OhFyYwHSDhl+5FKpOlUsIBEEpFAKruLgzgFIO7qvQeOrsypaQCIISCJkxQJwBgcjEbXt\nlkJAmqBUViyF3FiAtMO+BhWlxFIISBOUzYoNbQMApJGW9zUwddZRIhCQrqBUCSzMLQCkXdSTvDUl\nxrHlbEJAmsDCXhSMkADSy2yaoeXbitZUVboVU4CAxATkxgJIQTwYSemGKEB1C0gEAUnU6vYpPrDV\n8twCgBQmJsM1+ainopqqIgSkyxT8l6NZ5ciNBUgvLc89qKimqggBaYLiz1BIgQJIOy0fjKSulFgK\nAYkQZrJiMWUHIAUN3lnqqqkqQkAiROmsWErj24EAJKLNg5HUVVNVxG5ACgQCPT09Q0NDsS87derU\n4OBgiu/FQiSg4zMWWgLyiNvDE7wFIDZtHozEQh5LEhgNSE1NTWVlZWvWrLHb7Rs3bpzssnPnzj3x\nxBNNTU1peVNlnybmTRxCgX0NmhC3hzc0NJSVlb388ssOh6Ompkb+FnJDmwWE1LiARAjRK92AKPx+\n//r163fs2FFSUnLmzJklS5bcddddpaWlkVdWV1dfddVVqb8jfZpQ1sSpSN4AKVS6KSCxuD28p6dn\n48aNO3funDdvXl9f3w9+8IOysrKFCxcq2Gb1Eg9GqiIFSrdFPir9JGFxhNTR0ZGbm1tSUkIIycnJ\ncTgc+/fvj7ysqanJaDQWFxen632VHd5qeX+q1sTt4UePHi0tLZ03bx4hJD8///rrrz9x4oQybeWC\n1o53UWNKLMViQDp9+nRe3uXBZm5u7sDAQNg1fX199fX1L774YlreUfCOqG6yFdQrbg9funTp1q1b\n6Z8FQXC73UVFRbI2kS9aOxhJdTVVRSwGpPHx8YyMDPFLvV4/NjYWesHFixdfeeWVdevWGQxxoojF\nYrFYLFu2bIl9meALKP6PRyOi1vYCacTQ0FD3JT6fL24PF3V1dT355JMrVqy44YYbol6QYA/XOK0V\nL1ZdTVURi2tIWVlZgcDlZ5mRkZGsrKzQC+rr641G43e+853jx4/7fL7s7Oy+vr78/PzIlzp27Fgi\n78jIfKtNhUNsSMShQ4feeOMN+udVq1bF7eGEkLGxsTfffPOPf/zjr371q/Ly8sleOcEernHivgbF\nnzvlobqaqiIWA9LcuXMFQRC/9Hg88+fPD73A7/d7PJ41a9YQQk6dOpWZmTk2NrZu3bpU3lTxrFhK\nUzPd2uFwOBwOh/hlW1tb7B5OCHn++ef1ev3u3btnzZolTyO5p5F9DepdQCJsTtlZrdbR0dHGxkZC\nSHd3d1tbm91uJ4Q0NDR0dXURQlavXr3rkocffvipp55KJRrR4a2yWbEUcmM1Im4P//jjj3t7e2tq\nahCN0kJTByOpsaaqiMWAlJ2dvWHDhk2bNpWWljqdzpUrV9IV3c2bN0fdbpcidmIAcmM1Im4PP3Dg\nAN3IcMslzc3NSrda3cwmZc7elJ8aa6qKWJyyI4Q4HI6DBw8ODg4ajUa9fqKR7e3tkVe6XK60vCML\nDxRibqxKOxMkLnYPr6qqqqqqUq51HDIbDbSAEAt3uqTqOvtVuoBE2BwhUTqdLicnR7xXpcNO6g/i\nkKbI1sOBaOZgJLoAodIFJMJyQJINOyNcbdY4AZCBRhLPVVpTVYSAhKxYAP5p5GAkldZUFSEgEcEX\nYGTPN3JjASTF/UY7ldZUFSEgEcEbYGcHDnJjASSihYOR1L5IhoBECDNZsRT3D3EAiuD+YCRVp8RS\nWg9I7GTFUsiNBZAI95uG1FtTVaT1gMTapz9yYwEkwv0arXprqoq0HpAodp4pcG4sgHT4LiCk3pqq\nIq0HJNbyElT9dAPAOLqvQelWSIKDBSSCgMROVizF/TQ3gILovgYuZ+1UXVNVpPmAhKxYAM3g+IGP\ntWfr5Gg+IDGTFUtxv+4KoDgu7y9V11QVaT4gsZQVS6n9GQeAWbwejKT2mqoirQckBpmNBv5uGABG\nsPYAmhZqr6kq0nRAov+KrD1WcHnDADBCPBhJ6Yakk9prqoo0HZDY7JT0huF1cyqAsrg8GEntNVVF\nmg5IFGvjXOTGAkiHy4ORuImvmg5IbHZKDsbdAMwymwxmk4HN2ZHk8JESS2k6ILG8c5+bRx4A1nC2\nb4iDmqoiTQckNnE5pQDADs4ORmp1++iwT+mGpIGmA1Kr28fggXjmiTUkBCQASXB2MJLgDfAxPCIa\nD0iE1T3WZpNB8GJTA4AkeCogxNMCEmE5IAUCgZ6enqGhockuOH/+/PHjx8+dO5f0WzC7TsPZHDdE\nFbeHU6dOnRocHJSnSRrBU4EuPmqqihgNSE1NTWVlZWvWrLHb7Rs3boy84Ne//rXNZnvxxRfvvvvu\nHTt2JPEWbGbFUmyO2yCN4vZw6ty5c0888URTU5OcbdMCbgoIsbwzKwl6pRsQhd/vX79+/Y4dO0pK\nSs6cObNkyZK77rqrtLRUvMDtdm/fvv1//ud/5syZ89lnnzmdzkceeWTOnDlTeheWn49obqzSrQCp\nxO3hourq6quuukr+FnLPVmh0NXuUbkUa8FFTVcTiCKmjoyM3N7ekpIQQkpOT43A49u/fH3rB2bNn\nf/azn9EIZLFYMjIyLly4kNx7sfxkwXLIhFTE7eFUU1OT0WgsLi6WvYH84+NgJG5qqopYDEinT5/O\ny7sc83NzcwcGBkIvWLhw4YoVK7799tsPPvjgmWeeefLJJ/Pz86f6LnTulS5vsoazXUAQJm4PJ4T0\n9fXV19e/+OKL8jZNK/jY18BNTVURiwFpfHw8IyND/FKv14+NjUVeNjIy8uc///ncuXNnz549f/58\n1JeyWCwWi2XLli2RP2J/7hWzdtwYGhrqvsTn88Xt4RcvXnzllVfWrVtnMMTpnzF6OMSl9hESNzVV\nRSyuIWVlZQUClz+LR0ZGsrKyIi/LyclxuVwXL1587LHHmpqanE5n5DXHjh2TsKGSQW4sZw4dOvTG\nG2/QP69atSpuD6+vrzcajd/5zneOHz/u8/mys7P7+vqiTgOotIcrjo+DkbipqSpiMSDNnTtXEATx\nS4/HM3/+/NALtm/ffubMmV/96leEkGnTpt1yyy1ut3uq78JmViyF3FjOOBwOh8MhftnW1ha7h/v9\nfo/Hs2bNGkLIqVOnMjMzx8bG1q1bJ1d7QR0Eb4AUKt2ItGJxys5qtY6OjjY2NhJCuru729ra7HY7\nIaShoaGrq4sQMn/+/J07d3Z3dxNCTp48uXfv3ttvvz2JN2J5dzVyYzkWt4evXr161yUPP/zwU089\nhWiUdmovIMRZSizFYkDKzs7esGHDpk2bSktLnU7nypUri4qKCCGbN2+mm5Huu+++v//7v3c6nTab\n7Yc//OGPf/zjxYsXT/VdGF+hQW4sx+L2cJCB2g9G4qmmqojFKTtCiMPhOHjw4ODgoNFo1OsnGtne\n3i5esHr16l/84heDg4Pf/e53Q9eHEyRMbLHjZzEQ1CVuDxe5XC5ZW6YZal+p5ammqojFERKl0+ly\ncnLEezVSRkbGnDlzkohG5NLxd/MY/rek8wlKtwIkFLeHg6TUfjASTzVVRewGJElNjJAYDkiUeu8W\nAPapd2KcywUkotmAxHJWLIXcWACp0XkINU5FcFZTVaTRgMR+ViylxlsFQC0m9jX41LedVS2fYFOl\n0YDEPrWvuAKwT70FhDirqSrSaEBiOSuWQm4sgNRUejASfzVVRRoNSITtrFgKubEAUlNjASH+aqqK\nNBqQVLE2o94tQABqocb8Cv5qqoq0GJCQFQsAlBoPRuKvpqpIkwGJ+axYSo3PbgDqosZ9DYI3wOvz\ntCYDkkqyYinEJACpqWiERFNi2X+eTo4WAxLFclYspd4kCQC1UN3BSFzWVBVpMSCpa0kQIyQAENGU\nFbV8fE2VFgOSWiA3FkAG6joYSfAG2E9ZSZoWAxIt2650KwCACSqqG8lrTVWRFgMSUckMrNrL4wOo\ngoo22vFaU1WkxYCkolUZXjd3ArBDRU9+tJEcT/BoLiCpLitWRft/AFRKLVVRWt1+XlNiKc2dVqmW\nrFjKbJrR6lbBTAKASgneQH1nP73LCl7/lH6z4va8qsUFirYrCo5rqoo0OkJSCzqSU1ebAdSi1e2z\nv3PY1eyhww6z0WArNNoKja5mT8Hrn7I2icdxTVWR5gISpZZ/1HkTh1AgNxYgzao/8ti3fUYIaam8\nraq8gBBiKzTWOhfUOhd41pUSQpY3Hqn+yKNwK0OoK4EyOZoLSG2qmgGbOBUJIySAtKr+yONq9tgK\nZ3vWldoKjWEHI5lNhpZnF1bcnudq9ixvPKJoSy/juKaqSHMBiVLLUwZyYwGkQKNRS+VC8TthBYTM\nJkPV4gJb4Wx2Ju44rqkq0lxAErwjaolGACAFml667MrRRtTi+lWLCwRvgF6vLL5rqorYDUiBQKCn\np2doaCjpC6ISfAG1LCAR1R6xDIlIpAP39/cPDw/L1iSNqO/sJ4SETX9FPRiJTlHUMxCQ+K6pKmI0\nIDU1NZWVla1Zs8Zut2/cuDHygoaGhrKyspdfftnhcNTU1CT+yqpbj7FxvctTs+L28BMnTjz00EPL\nli174IEH1q9fL38LORYjmyesXgM7tcD5rqkqYjEg+f3+9evXv/3223/4wx/27Nnz/vvvf/rpp6EX\n9PT0bNy48f333//ggw927dr1m9/85vDhw4m/vurmYVm4HyCN4vbwYDD485//vLKysrm5ee/eve3t\n7YcOHVKqtZwRvAGzyRCZzUMHQ5GzEXRmT/EHWb5rqopYDEgdHR25ubklJSWEkJycHIfDsX///tAL\njh49WlpaOm/ePEJIfn7+9ddff+LEiURemfY2dc3Dmk0zFL8ZIL3i9vDDhw/rdLqlS5cGg0GDwdDc\n3FxcXKxQYzlEk2HDvmk2GVzlBa1uf9iKUZvbT2OYjA0Mx31NVRGLAen06dN5eZcH1Lm5uQMDA6EX\nLF26dOvWrfTPgiC43e6ioqKoL2WxWCwWy5YtWyYuVuEnO3Jj+RO3hx8/ftxisbhcrpKSkjvvvFPs\nwJHCejjEFWMWbpk1z2wyVDdfkXtU19mv+GZrjSwgETZLB42Pj2dkZIhf6vX6sbGxqFd2dXWtXr16\nxYoVN9xwQ9QLjh07FvlNdf27irmx3E8fc2xoaOjkyZP0z/n5+XF7uN/v37dv38qVKz/55JOvvvrq\n6aefNpvNjzzySOQrR+3hENsya97yxiOtbl/YR4HZZKgqL1jeeGR545Fa5wLCzNBE8I4Q9SSrpILF\ngJSVlRUIXH5+GRkZycrKCrtmbGzszTff/OMf//irX/2qvLw8wVdWY0KPFnoh9w4dOvTGG2/QP69a\ntSpuDzcYDFdfffUzzzyj0+kWLFjw2GOPtbS0RA1IkAQah6o/8tgqw59NK6x59Z39rW4fDVeMDE24\nr6kqYjEgzZ07VxAE8UuPxzN//vywa55//nm9Xr979+5Zs2Yl/sp0nK6uj3jxsBbF7wpImsPhcDgc\n4pdtbW2xe7jZbM7MzNTpdPTLzMzMCxcuyNJSTaDLRbQKAx0Jhap13lTw+qfVH3mqiafV7XeVF2AB\nSTYsriFZrdbR0dHGxkZCSHd3d1tbm91uJ4Q0NDR0dXURQj7++OPe3t6ampopRSOCrFhgQ9wevmjR\nouHh4Q8//JAQ4vV6d+7cabPZFG0yb5ZZ81zlBXWd/ZFFVM0mQ4U1r9Xtp9GIkbLfGnke1QWDQaXb\nEMW+fft++ctfTps27dy5c88999yKFSsIIXfeeefjjz++evXq6urqxsbG0Fn4t956K3LizmKxhM2w\nF7z+qdloCK0Xogq6F/eFlTkBtYvdwwkhn3322T/90z+Njo4ODw8/+uija9euDe3wVGQPhykRvAH7\nO4dJxMc9HZSYTQZaZVVZyxuP1HX2B//NEf9S9WM0IBFCgsHg4OCg0WjU65OcV4y8XXUv7quw5kUO\n0hln33aYEIKAxJlEerjf7585c2ZkKKIQkFIneAOtbl/YzrqK2/Pu/d5s+7bPWPi4KHj9U1qGXNlm\nyIPFNSRKp9Pl5OSk/WVVlxVLITeWP4n08NmzNbFyoCCzyVBhyou6ZYBWVo3cjCczwRsw367KT60k\nsLiGJBE1ZsVSyI0FkF+t8ybBG1D2+AmN1FQVaSggqfczHbmxAPIzmwy1zgXKVvum+841suebaCog\nUWrcrIJzYwEUYSs02gpnh60wyYnWVFXq3eWnoYCkxqxYClvVARRBj+lTcOJOIzVVRRoKSGrMiqXE\n3FilGwKgObZCY4U1r66zX/5jyTSVEktpKSAhKxYApq6qvIAQIv8giZHCRXLSUkDyqfVEepwbC6Ag\npXY3aKemqkhLAUnNs7Ga6pQArBF3N8i52VU7NVVFGgpIqmY2GpAbC6AUcXeDbDvuNLiARLQTkOhz\njXr/ddU7tgPggyK7GzS1gES0E5DUvgBjNhoEbwC5sQAKknN3A93RoLW5eq0EJErtjxvIjQVQkLi7\nofojySfuWt0+rS0gEe0EJPVmxVL3fk+tk40APKG7G+oO9Us9XSF41borOBVaCUjqzYoNhSk7AGWZ\nTQZadFXS3Q1aq6kq0kpAUjtarEHt4zwADtAjZSXd3aC1mqoirQQktdcoNE/UV0VAAlCe1Lsb1P55\nlTStBCSi/p3TZpOBZm4DgLKk3t2g6iz+VGglIHGw+oLcWAB2SLe7QZspsZQmApLas2IpbT4xAbBJ\n3N2wvPGL9L6yBmuqijQRkNSeFUvR3FilWwEAE+juhla3P72fMPTV1L4lODmaCEgUH//AfARXAD5U\nlReYTYb07m4QvAEN7q+jNBGQ6G5punNavZAbC8Aas8lQVV6Qxt0NWl5AIiwHpEAg0NPTMzQ0FOOa\nb7/99uTJk3Ffio+sWAqzdtyI28OHh4ePHz8e+xYAxVVY89K+u0GbC0iE2YDU1NRUVla2Zs0au92+\ncePGyS57++23t27dKmfDFITcWJ7E7eHvv/++zWZ7+eWXy8rKXn/9dflbCIlL4+4GbdZUFemVbkAU\nfr9//fr1O3bsKCkpOXPmzJIlS+66667S0tLQa2pqatrb2w8fPvzDH/4w7gvykWWG3FhuxO3hw8PD\n//zP//zb3/72tttuO3Xq1P333//www/fdtttCrYZYjCbDK7yAlezp9XtS3Fwo82aqiIWR0gdHR25\nubklJSWEkJycHIfDsX///rBrSktLn3vuuR/96EcJviYfe6aRG8uHuD1cp9Pp9fo5c+YQQmbNmjV9\n+vTMzExl2gqJWWbNS8vuBm3WVBWxGJBOnz6dl3f5GSE3N3dgYCDsGqvVevfdd8+bNy/2S1ksFovF\nws26C3Jj+RC3h8+YMWPt2rWVlZXbtm176qmnli5detNNN0V9KdrDt2zZIm2LIR5xd0MqMUmzNVVF\nLE7ZjY+PZ2RkiF/q9fqxsbHkXurYsWOCN1Dw+qdafugAxQ0NDYm7b/Lz8+P28AsXLnz99ddnz57t\n6+sbHR3t7e31er3XXHNN5CsfO3ZM0pZD4iqsefWd/a1uX9ITd5qtqSpicYSUlZUVCFweB4yMjGRl\nZSX9avRQOz4eOmyFRm5Ge5py6NChf7zkwIEDcXt4V1fX7t27d+3a9dprr/33f//3+Pj4e++9J3ur\nYcomTqZIdgs4H6vdqWBxhDR37lxBEMQvPR7P/Pnzk341+gnO066V1BdOQWYOh8PhcIhftrW1xe7h\nPT09+fn5M2fOJITodLri4mK32y1XYyF54u6Gus7+JAY6gjeg8VubxRGS1WodHR1tbGwkhHR3d7e1\ntdntdkJIQ0NDV1fXVF+Nj6xYCrmxfIjbw4uKiv7yl78cOHCAEOL1ej/88MPi4mJl2wwJorsbkji+\nT+MpsRSLASk7O3vDhg2bNm0qLS11Op0rV64sKioihGzevDlyu11cPGXFUpi1U7u4PfyWW25Zt27d\nqlWr7r777vvvv3/RokVOp1PpVkNCkt7doOWaqiJdMBhUug3RBYPBwcFBo9Go1yc5r2ixWI4dO7a8\n8UhdZ3/w3xzx/wfm0Q0arvKCqsUFSrcFUhW3h1+8eJFeMH369KgX0B4uZRshSfZthwVfoNa5IPEA\nU/D6p4I3wMcnVdJYHCFROp0uJycn6Wgk4mmdELmxPInbw6dNm3bttddOFo2AZUnsbtByTVURuwEp\njfjIiqWQGwvAPrq7odXtpytDcWEBieI/IHG24oLcWABVSGJ3g8YXkAj3AWlizzeyYgFAXmaToda5\nIMHdDRqvqSriPSBxlBVLITcWQC1shUZb4WxauyH2lRqvqSriPSBxlxVLISYBqMKlkyniDJI0XlNV\nxHlAovjIiqVobiwd+QEA48SJuxi7GyZ2NCDtnfuAxOvMLEZIAGpBJ+5i7G5ASqyI84DEH5wbC6Au\nZpOhanGs2g085UqmiPOA1Or28Tc8AgB1sRUaK6x5dZ39UXc3CN4AT7mSqeA8IBHuBsJmk8FsMsTd\ntAMATKkqLyCERA6SkBIbisXjJ9JI8AZIodKNSDfsxgFQHbq7YXnjEVpPqO7QxB4Huh7c5vZj2zfh\nOyCNzbiGcPrxjWINAKpDdze4mj0k5FjYOm8/IYTO5k2pGCuXeJ6yG8u+hvCVFUuZTTOwyw5Adeo7\n+1svHVJe61xA/6NfetaVEkLE8ZNm8RyQxmdco3QTJEHHfIhJACpS/ZHH1eyxFc4O3d0gLiCZTYaW\nZxdW3J7navZM9SAlnvAckCj+hsDzJg6hQG4sgGrQaNRSuTBydwP9jKK7w2mpIcVaqTSeA9LwNfOV\nboIksJEdQF3oSGiZNY+E1G6o/sgTmblPM5YSPLSCPzwHJIq/j2+aG9v2pV/phgBAQuo7+0nIRga6\nu6HuUH9kTVV6d9cjIPFnLPsa/qIRAKhO65W7us0mAy26GllT1Wwy2Apna3YbLc8BaXzGNVzu+aZR\nVssTzQAqIngDZpMhLPXVbDJUWPPMJkNkTVU6s6fNXUtc5yFlX8NrQQ5UvgJQEcEbqO/sD5udo3u+\nI7W5/TSGydI0tvA8QiKcZsVSmh3UA6jLVGfh6iJCl3ZwG5DojBZ/WbEUcmMBVGSZNU/wBhKZZtd4\naTsVT9kFAoGTJ0/OmTNn1qxZYT8SBKH+g4OEzDkhCET6Zw1BEFpbW0+cODFv3ryKigqp344QQoZO\nEUKq/337vbd932azSf1ugiDU19fPmzePECLPX7Curo7+Pm02m9lslvrtxL+g2WyW4feZuG+//dbn\n8/3d3/2d0g2BlNBMo+qPPLbKOGmRWj8bKahOu3btuv3223/wgx8sXLhww4YN4vc9Hs/EB8rNi8kL\ne8msXJvN5vF4pGuJy+UK/X2azWaXyyXd23k8HrPZPPG3u64oibfbvHnzlK6P/AtK+vtsaWmR+fcZ\nFoGk7jBT8sYbb7z88suT/fTGG2+UszEqMtVOLgPXnq/IC3sr3v1isgsKSVE1AAAGoUlEQVQ834zY\ntnaRF/a69nwlXTMY/M2EUmVA8vl8t956a3t7ezAYHBgYuOOOOw4cOEB/dPnDZdGTNCDRjxiJWhL2\n6Sl+hra0tEj0jhPDheuKaECi7zilj+wpfYqFRSPxLzj1hifE44leyEu632fU8ZB0HSZxmzZtcjqd\nN954IwJSEhj8zXi+GaExyfzagZYvvWE/bfnSa37tgNTRKMjkbyaUKteQOjo6cnNzS0pKCCE5OTkO\nh2P//v2EkNbW1tbW1omLZuUSMjG1dcX302r58uWR3xQEobq6Woq3q6urEwSBEEKGThNCyN8Vid+X\n4u0IIVEDkiAIEr3jZL+3qL/n1E3WMaTrMIkrLS197rnnfvSjHynbDEgXWhlILKIa9p9922eEkJbK\n26oWFyjdUiWpcg3p9OnTeXmXV4Zyc3NPnDhBCJn4sKZm5dJoRLW1tUmxNnDFOybw/RS1tbVd8fXN\ni8l1/4cQIhBi33Y4wRc5uWh1ghcLwgny2FtRf1T9xdX1Cb9j4lrP3xr1HYWp/AUTJwi+K97u/74g\n/lGiDpM4q9VKCPn8888l6kugCFpEtdXtq26+YjLAVV6g8VBEqTIgjY+PZ2RkiF/q9fqxsTFCCA1L\nE0I+XAgh27dvb2hoSG8z6JtG9de//tVisaT37Qghp05dCrFDp8jBehqNqAMHDuinT0/wdf5fR0eK\nLfnrX3tPD5xO8UUiTZ+uHxsbj/qj1NscaXzyf0G6g0M2Q0NDJ0+epH/Oz883GhNa077jjjuk6GZ8\nYPw3k3nllw3/RdL88TSJO+64Q5b3SZIqA1JWVlYgcHnT88jISFZWFiHk3nvvnex/effdd6V44LXb\n7VHndn7yk5/U1tam/e3q6uouT14d/C9C/kv80XG62SGtBEEoKIj+1LajtlaK7XaT/T7pToq0v11r\na6vdbo/6Ixm29oU6dOjQG2+8Qf+8atWqJUuWJPJ//eY3v5GyUQCyU3oRKxmtra2lpaXil5WVlf/+\n7/9O/yzzGnXUTQ2EEIm2aXkmiTrS7UPDpgZlvfPOOzE2NQBwRpUB6fz588XFxe+++24wGPz8889v\nvvnm//3f/6U/8ng8YZ+hUn+4RO6Klu7TMxgtJkm9K1rmv2BLS0vYX7C2tla6t5O/w0wJAhJoii4Y\nDEZ9JmXcvn37fvnLX06bNu3cuXPPPffcihUrQn8qZqouW7ZMnrRKORNj6dsRQuT/C957773y5OGK\nf8Gqqiqp347I3mESt337dkEQ/vVf/1XphgDIQa0BiRASDAYHBweNRqNer8qVMAAACKXigAQAADxR\nZWIsAADwh8+AFAgEenp6hoaGlG6IwmL8Hv72t7+dC4GB8rfffitmArEPPVyETp449js5h6svTU1N\nr776an5+fm9vr9PpfOmll5RukTJi/x5qamree+89g2HieI7f//73OTk5SjSTFW+//bbf71fF9gH0\ncBE6+ZSooJMrucVPAjHqrmpK3N/Dz372s08++USh1rElkTKm7EAPF6GTJ04tnZy3KbvJ6q5qTdzf\nw9GjR7///e/39vaG1rzQJnWVMUUPF6GTJ04tnZy3KbvJ6q5qTezfwzfffHP27NmKiorh4eEzZ878\nwz/8w/PPP69EM5mgrjKm6OEidPLEqaWT8xaQJqu7qjWxfw9nz5598MEH165de+2113Z3d//kJz8p\nLi4uLS1VoqUwNejhInRy/vA2ZTdZ3VWtif17uOGGG958881rr72WEHLzzTeXl5e3t7cr0EqYOvRw\nETo5f3gLSHPnzg0dk3o8nuuuu0655igm9u+ho6OjublZ/PLixYsXLlyQs3mQNPRwETo5f3gLSFar\ndXR0tLGxkRDS3d3d1tY22fkCfJvs99DQ0NDV1RUIBKqrq3t7ewkhX3755b59++677z6FWwyJQQ8X\noZNzSOltfum3d+/ekpKSRYsW3XLLLe+8847SzVFM1N9DSUnJW2+9FQwGa2pqFi5ceP/99y9atIjW\nTdc4FdXVRg8XoZNPCfudnM9adkHUXSWExPs9jI+P+3w+jacKqhR6uAidnCd8BiQAAFAd3taQAABA\npRCQAACACQhIAADABAQkAABgAgISAAAwAQEJAACYgIAEAABMQEACAAAmICABAAATEJAAAIAJCEgA\nAMAEBCQAAGACAhIAADABAQkAAJiAgAQAAExAQAIAACYgIAEAABMQkAAAgAkISAAAwAQEJAAAYAIC\nEgAAMAEBCQAAmICABAAATEBAAgAAJiAgAQAAExCQAACACf8fblhnwyo5nXgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = 0:1/5:1;\n",
    "u = zeros(length(x),1);\n",
    "u(2) = 1;\n",
    "figure;\n",
    "subplot(1,2,1); hold on; \n",
    "plot(x,0,'k.','MarkerSize',18); \n",
    "plot(x,u,'-','linewidth',1.2);\n",
    "subplot(1,2,2); hold on;\n",
    "u = sin(2*pi*x);\n",
    "plot(x,u,'-o','linewidth',1.2,'MarkerSize',10);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAGkCAIAAACgjIjwAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAA\nB3RJTUUH4gkHFSsYeG5ORQAAACR0RVh0U29mdHdhcmUATUFUTEFCLCBUaGUgTWF0aFdvcmtzLCBJ\nbmMuPFjdGAAAACJ0RVh0Q3JlYXRpb24gVGltZQAwNy1TZXAtMjAxOCAxNDo0MzoyNMJBO+8AACAA\nSURBVHic7d1Nkus6cobhpMNTb8CT6m14pLOn6znFBfSepFFvw5pcL8ALkAdgQSAIgCAIUEnyfaKj\no65Kf6Uj8VMCSaB7v98CAMC3/du3nwAAACIEEgBACQIJAKACgQQAUIFAAgCoQCABAFQgkAAAKhBI\nAAAVCCQAgAoEEgBABQIJAKACgQQAUIFAAgCoQCABAFQgkAAAKhBIAAAVCCQAgAoEEgBABQIJAKAC\ngQQAUIFAAgCoQCABAFQgkAAAKhBIAAAVCCQAgAoEEgBABQIJAKACgQQAUIFAAgCoQCABAFQgkAAA\nKhBIAAAVCCQAgAoEEgBABQIJAKACgQQAUIFAAgCoQCABAFQgkAAAKhBIAAAVCCQAgAoEEgBABQIJ\nAKACgQQAUIFAAgCoQCABAFQgkAAAKhBIAAAVCCQAgAoEEgBABQIJAKACgQQAUIFAAgCoQCABAFQg\nkAAAKhBIAAAVCCQAgAoEEgBABQIJAKACgQQAUIFAAgCoQCABAFQgkAAAKhBIAAAVCCQAgAoEEgBA\nBQIJAKACgQQAUIFAAgCoQCABAFQgkAAAKhBIAAAVCCQAgAoEEgBABQIJAKACgQQAUIFAAgCoQCAB\nAFQgkAAAKhBIAAAVCCQAgAoEEgBABQIJAKACgQQAUIFAAgCoQCABAFQgkAAAKhBIAAAVCCQAgAoE\nEgBABQIJAKACgQQAUIFAAgCoQCABAFQgkAAAKhBIAAAVCCQAgAoEEgBABQIJAKACgQQAUIFAAgCo\nQCABAFQgkAAAKhBIAAAVCCQAgAoEEgBABQIJAKACgQQAUIFAAgCoQCABAFQgkAAAKhBIAAAV/v3b\nTwAANum6Lvar9/u95zPBRlRIAI6ns4bu9r6ZC90fzP+6oZv73rPGgo5vEAAOwc0Smz3P7rl4iev5\nfMpD5D7+JwdAVQgkAKrZHJoHzLN7mgvtD+6vxp/u4bu99Td7NQ6DSjCHBECdTzF0Dxc6YrLkIc/h\nKXeRuzyHp/vb2/tmUmpMprvIbwhN3MeHI5M04J8BgAqTEJomxzjO5rqPV5vXRuNNBieHfmMpfU0O\nhl9HhQTgayYtBg+xqeOXO/1NfqPk2T1jZZN7qzHS+vEmdmRPIrF0e9+ok76OfwAAe3OLISswnjbj\n1UMmnHLuwb1hYMJpeJqbmwFAjorfQiAB2FXXdbf3zauBjHQmfQod2yk3G9zLubnIOPPk/qdJI9vp\nwIHxKwgkAPsxtVG0TyGUUvJH5I+IfMopGxux+4nxhuzcWCKTNCCQAOzHlEerbmLH5dxiqCCNYre1\nbXhjFDmZJHSE74tAArCT4jQaZ3dEZEN55N1t8Lwlt1qyl3OQ3A1LBwHYw9o1e57d040NkxAmmTam\nkfw2gn9OnjVRdBfvQvntvtvyWMhHIAFoLj11NGd7u91iZfR7eXjCKZtZ7G6SSf1tnFK6i5dVLIK3\nD85DAqCLO0znplFg4O63Xbv4seyCDt6FMh3WMz9wolJrvL4A2lpVHsXSyP42sKLd5ljyOsI/E0iD\nfxIuU0pNUSEBaCg/jdwlGJ7PQBqNV5sFldcdVyC2ttCtvz3l6f6KBR2aYg4JQCv58y62MBr/+xG+\nTiLYxn6H0oklO6XkDgyKM7Fk+x3sNktlD4QEoh5AE4ltIzxeGi0O1qWLoY0jeH40zrrAxZlk4vhZ\nF0N2AJq5R7cjMuYDZVtG3oyN/Q5jMXSfdJzbZ+Ut0srwXV28mgDqyymP0rWId835cqhZi7GWxlKg\nnWH2iG5rOAfSKphDAlBZWRq1cOtvZXNLZurIO0tpvgeg2+lQ5QlfHBUSgJpsW12iByGYRvnlUfr6\nMfnVkrdYw0QoRMdda6mTNiOQANSULo9ihdHaNErcJG3eoTDnLQc+uSR0E3diiSPqFgQSgGrS5dHa\nNJKlVu+cTEqUO1ntf49x84tEj9/8ITiulqHLDkAdZWmUsGoR1cVyx7/+8MyqsR6fk5DGGaP5mkah\nhcNRgAoJQAXjPrCh9Q6CF35+u6Y82lLuxB5dIhNL7ga1t5v/py2cC8UpSkWokADUYY7gXmakC6N0\nGsksgYIL2W1aXDVy0tIkCx8it/HRE3XSxJ2VWEvwkgHYKjZYtzaN5gVQ1iJ4m8+ltfcjkQ0Ag0VS\n+qEpkgpQIQHYZBysm53os5BGv8un+rvk/d52vgdEa59q6T4LwkegSJJ0nUSRtB6vF4BydurIPVIv\n9i8EQ8j77TycEmoVSfL75GU2sWSLpEAh+HzebuFZKA6wq1AhAdjgLvlp9MmhROU0m4iK9Uq0E55Y\neoxFkswHJ51fTVAkrcSLBaCQu9fRJ0jufneA/Xm85lJ/mkSyZ/s5SYuCi+bJdMul4HBi4sRejrH5\nqJAAlHD7vO3/u2sWGKsWRQ326bl3teeU0vigTrWUeD6xySRWBF+FxVUBrOZOHbmX2xNI7f/y7zMn\nbNyN8vxfzRY/XSvxBMwirWICNRKZ258AqJAArOOlUf4ET/ESQa5G3XfpRfPMD33fi8gwDPbPz2m6\no0jKRyABKFcljbw7zLm3YH22cLLqGm6tY3LI/c9hGGLP5NbfAk139+3P6BIIJAAreH3eUiONysqd\nit13458TDyFP3/eDDJ8JpKWmu1tPkZSF1whArkkazbZmSMgfrCsIJ/8e1hRJYwjdRe4LITQ3dIPb\nUphuujMvF8fbNAIJQBaTRuIexEUkZxOH7PKoeGYopwPbe0rmh77vh27o3+uiyFqXSbSAL2HIDsAK\n7mmwdlGD1PUbDNbN5XRgB0fktqSRvathGAIN3/OnwXmySwgkAMvmSwSN8zfJWmSfNDKCU0r500Jl\n+vdYYJk2h3TT3a2/Pe/0hacQSAAWdF0nj13XOS10F/ldu8hEUTqEtpdH1rwj3FzuN91RJCURSABS\nuq6Tu9xut6d8lmNIL6lgrC2PMgPv+XzKI/yr8eF6eXbPWkmzyBZJYofv4k13FElpZDWAlG7o7AJu\n+SfDFgzWzZvlgrfNaVgwrQoikoiliuWRd2/DMMgfkT/h7WVpt0ugQgIQZcqjcRBsuqhoojxa0VY3\nLXfcENp4fquJh1jq1E0j/6HN8N19CPZZ3Pqx1sQcFRKAsK7r3MO6aXEWGc/aMf8LNrOlO9zG24qI\nhBumtzDlkXtJMHtaBNL8PoduCJ4gRZEUQ4UEIOIu4qSRV3MMMniL6Fjz0bbqwRM0TyP5neMRZ/iu\naXk0f+h5050pkuhumCOQAASY8shLI8udxvduGDwpx9g/jYz08F0twRfKVpaf8brI9rIQtp8AMGcH\n68Sbrt+rtmjBZMP+f0L/7m0buoiYObNbfzMt4Hs+E/0IJABhi41qQX3fBxvkUvs7bK6cEuXRnPm7\nWrDDg/PLbSaxbVICgQRgwh2scy+fTyPl+9ZgnctGhS2VGj2fIJtJ4z7opsGBImmKQAIwdR//v3h0\nK1YktZCTRu5Inc3aRrEUK5LG3/a9eXGok4IIJAAfXdeNZ5XeJ4N1W6Ze0uXRluIpM41M/NhL3L67\ndH4Uuocv/jzubyaZK1MkWQQSgJEdrJP41FFmMtki6bsr4OV0ZFQvlew6EbHHEvv6/Kn1mCdBIAH4\n2O000u0Wy6PE0/YKo+ojeInCy82k/taLUCR9EEgARGbrMmy/w77v25VH6TQqa+9uN7E0f6BPh0VP\nJn0QSABEZOxikNlgXXE+fZYaiitLrMU08iaNghIt2lViaaG7YZ5JYKUGACKfr+fp43h+OH16xOML\nN7RQazVVWyzK0mtSDVslUSEBcNdlCLhvuvPqLeCx8qhgmG6xxW5jtUSRtBaBBCBaCgzdUHCszAyG\ngvG6RBrlDNOVadUdPs+ky88kEUjApS0M1t0/P2bGTGCB0UpFUjqNVj2lz3PLTpqyUimnCNt5zQjN\nCCTg2u4ikTQqK492tr0ncFUmtWjDcxvBL14kEUjAdZkNYRvt5O3aXiQFd977yglSa2MpJ/CokwwC\nCbio9DfxgmP9qpusmkCK7QNbK40K8mBVLK0owq5dJBFIwKWtbeMuf6DSIilzV/KYplVUxUG8SWjd\nN97ZURFIwBWZ7+AbF6zbeJMC7Ybptgya5cRS/sCdSd9rFkkEEnBRO0wd+Y+4fsVVtzxq2tu9XZXu\n8M89XHLgjkACLmft7FE6b9rVRvM0WnsPq25VpbMgUSpl3r/Zyu+aWDoIuJb0YF1rfd8PQ9ZB36bR\nruv31LB92aExuq63mBAVEnA5FcudksLlLovjdW4a7TlMV7H9OjixRHt3GoEEXMjCsnUrFQ/WZbbb\nbRwMLLt53cwobsMbB+4uNpPEkB1wFYuHNiW785nySOeugGW8QbzcHpDrlVNUSIfUdd2lvjehFg3l\nkSTPSXoOT7Mz0xfTqN1SqjaN8ld5kCu1gFMhHY99d15twhNN7XO2afqGJo1OUxgF2TTK7Xq4X+g8\nWQLpsO7s6IUVymaP2lUqpt3O37uvXhptfOb5A2vFd24+ucHqx33cvu8HGS7ySSeQDqbrutv79uye\nt/72vD9v79tF3qlo6uuL8Ty7pxyqt3sLm0YiEvzwBlPqCp905pCOxKzNXLCtGVDLljTybusu3KBw\npK7RTJKbRjHvmepPQycC6TBMGk2GOH7D6Tpznih3j/4mlTH37GuWGtOo6sZLZ+rQMxLje2dCIB2Y\nDScyCWld1xXuRN56g7676hn76kVSTnl0ZcwhqTaJGa88+r3QFElMJqFMfjFRt+w4yoJATbsbVnm/\n3+Z8jxN/zKmQ9DJtUeaT8PP+EZHn8/nsnvKQ53M8h8ONKOoktLP9oOzeg/J1uxuhPFpEhaTRfPnL\n1/D66X9E5CWvn9sYTuPvnDYH6iTMLTR838MXtysLWhcc1e+/SpFUJY1OXyQRSOp4h4+hG37eP6/h\n5V3NxJKIvB4vcfMJyJZ5nK11iD/KMB2+hSE7XWJpZMqjIPOrn9uP+Z88xoWFGL7DsvvyVaqkkY2i\ng6bRxu6GioN15263o0LSIrZLzeuZSqPRXV7dy8wz/dx+5C2v7tW/e/dde9YaH8vu4YsXkiZyq7V2\nXh5USQMCyhBIKgRH+T+Ddb+NCzZ1Fv28f7xPJuF0TRtXDCrep9X+7L7Z6m5+sbPimaTqvQwnnkki\nkL4vnUbL5ZHIT//zkpcXV24med9SbTid7w2NTDnH1oLN+swPsfeVOZIeN5MK0Fm3CoH0ZYk0Kri3\n13Pa+/AITCO7R5lu6OzIDB8bFIgVQzGHziTOSWqNQPqa2KSRFSuP/NSx/kfkH+Il2Ute7skfMltI\nWH7/i3A6n9ihP31UtRv2xK6zNoQ85mAqDdrtlKSF1bQ8si/jmRBI35H4kuiWR5Ps+SPykNfzZRu+\nJ24ybw0XZ+Bu8bNKOMGav1s2hpDHtoqpyo8ceook42RFEoH0Be7nMNCD9JDX8yUP8coje0rsWvMG\nh0WE04nlvBnm788W//QHHb7LzKQdZo/OVyQRSLuy7x77OZ9/DzWpY053Xevn/TPvxHs9XwWZZE3C\niVa9Y7mvvoW7mam0/1euOHynqnDZs5fhTEUSgbQf+2UwdmaGHazLbK7z/PT+gg5mxM8k3JZMsugj\nP5BYv8zie2Dnf80jDt/pGbg7WZFEIO3E/bx5PW/250Rn3aozkMyV57NNVTLJckd1CKdz+GKP8kGH\n74L2fxlPUyQRSM2lRyTcw/rYwvAQ+bP1QaO9DxvECjvvDErCSa3F5rrv/nsdK5MoklogkNrK/IB9\nBut+Fwpy++vyo2UctXuI/BEJHVsWi6ToQi/xHaa9bnKG9TQ40JHdVTylpCQb5Hu5fo4iiUBqaG0a\niYg8xCwU5IaQCSd3QijGXucl0Z4Ik0mx39Y9EBBOeigvj6wDTSnNi6RvvYynKZIIpCbKvuXFehlM\nwIzdd/HKyQzTue156VNod+5rIpzU0pNG1rGG7zQwr5i2f8e1Dv8HKLTqgzQ5DTa2NEOkoyGaNw+R\nPyJ/Ul0S5m63fODr7pFj8G6spRu6vp98c0+MuKp92TO/2H13vM5dMfK7r+QJAokKqbLqaWRF48eU\nUNPfpkftxuuUNt3V3WaNyqm609QWhxi+09bdcOgPDoFUzdphutgKqsHgiXXN2csnc07yMq0N+c3i\n+Zp+9ginFg5aHlmHGL47xCupH4FUx9oPjCkyPtnzEPkjr+dL/hFoWwguGpTocTAzSaZIsg8RvubK\nIqldGqX7LAinFg50DG23HmuZ+dtVySt59CKJQKogvVJq7FZu7WLyI3i1YJWTeZqRu5LQ6/WS//m9\n3LltZiZVGaZb21NunxiVUz7v3Rj7xz1QGhmx4btvfUnKvOZXHDeTCKRN5mvTTcSPs5M0WrlQUE4a\nmROS3Lv9+fmR3//ywinn5KSyjdo8q0qx2PUJpytrMXy3/e2qyqFbwAmkcu6/evF7154Jm3llWXOe\nbEwgnB6pr9Kxy8P3Hj+FNlN+/hFOy+6Byw5XHrnKMqniuXeHcNAiiUAqNB8YsT+n399+L8NjPBN2\nUdlqQMH1vydXsOH09v8K+5/BD3O7cZLCVckJp7zxukOnkeEWAZmjZ+1SR+HredwiiUAqMf+ClhlO\n/mDdUnn0mQFan0bzUbvUA81a+/IjtpaK3eTBcNJ21PgKhUfPMnZKyUugU5Y7BQ7a3UAgrZPT6rOi\ncnosl0fbh+kmrQ2Rk5nG+3+PZ8vOF0SxP2srjBbZu73aToPBXV9P9ocfoiMc+QikFQre+v5h/eEs\nN5dXvhREkZs69udYT7kvNJnUOpxan1c4H8Y54pdHKKQ549U+sQQCKdfGL2LzwTpZqn4SSzOkr/C5\nw9tvC9876yTZcWDwnUqIuuFUcZguv8lC80FkC/cteoXyyKBIOhMCadn2M/ICizI8xJZH8/VS3UW7\no7GUU+6s4U5TZZ6ctDGctixclH4y6due8ricdtY0wskQSAt2+PLlhtBktK17yaNy6gQFC7W1izis\nDafq2zItus5BWcmeCLuhSDoNAimlyrt8soKqsyfspPT5xyeWbDCY1VE3NjXYXrtY/3eLvWXT4eQO\n0+3WpHv6gzJH5P2d/k21PwIpbOMw3eQ4+/hkz5g6kYWC5tnwCaeiDWQXpdOoeDlwV+x8pqEbtp9C\nm+Nqw3RXK48MiqRzIJAC8nd6Df/COc4O3eBvoxdqrlssg+YbyM4vXzTv/15egmhNJi2WO+5dmUAa\n7jV3sgg+pSscjmMu/ufjcGh+9QVXbwxeM2dBUm+ILFiRhC/M2zkiM5xsCpq7XTtM527ll4rhPvqC\npLvpGp3ndKnDsXnfepF/nT/f2LNIuuDLuwMqpI/wSqmlw0rh7Y4e/pmwG6dwgpVTehSuIOrMyUnm\nx7K+uGqnEmc/olxpmI7RKpwDgTSyXzDtJdU/4d5gXa2VUq3PhNNspwnb2uA1ka/oKX8Xbnm+dhZq\nezhd9qvrxcsjYSbp+AgkEecLZq2v6uHyyNGit80KLuYt3qoN5pote8q3n/Ra8M9xzQOx58ovApl0\naFefQ8rsplsVTsE0csujxTQq23r89XyZ4Anq+34Yfv+KuxTcv8iKIqnpakDBf46rDdNZ3uzRldPI\n2CGQeJEbuXSFlP/Gzf+qnq6NqgzTvYbwOFvf98GlWodhMO0GtulguA/uneRvyJS/vWzTI0Lwn4MD\nhHCgFBGKpCO7biAVv2ULxpHs/E1+FCXKnUQzm8eURPPr++1Yw+evWAyn6tvLZgp3993Hv+XiR+HW\n3wCAfVxxyG772nRB4xHzITKrgV7DS/6EC6NwuVPjjFFbGIV/G9sHNi+cgmN32yeN1vaUX3aYzmW2\n1SCYPe2KJF7ndi5XITWt5f0lgoyHyB+Rh7wegewJH2fv0e2WM6XTSGZFUvD5JMJpXiet+pKeKHcy\n70E4LoiI8+2KVwMncK1Aaj2yHOiifjj/nzfatnH4JTZMVyAznCS+VXb4fqvUfxx/Hbwac8wkHdFV\nAqniMF3sOBvsrPOCwT2sV8kM/7ktFUbFAuF0n61N52q2Th3DdB7v/DnguC4xh1TwRWnVZEZwiSB5\nLEROLJyKK6SCNFqxTt0QHmRz7bNS6hXesZns2iK8JjHViyTegU2dv0JKvyNbbH8wFkahDuzJQ0QG\nxArUGqZLPI1gDCcW/WvUaMexYI7XJIGBu2M5eYU0XxBoIrkeaCZ/b/LZMN3ae3PLjqw5p4LCyA2e\n+6fQybyfnG666uFEGs113ck/v1XUDSTeh02d+Q1tBzSaLhkw6axbGqbLucNJ69rSnFMijTLLnbUj\nhAUjihvDiUmjGAIpU61MIo1aO/8b2saS1E6mSRptK4zc+4zu0TANJztMFwueFefPrtnxaONruDac\nOARgOwLpKM4/h+S+gRqFU600SvBSx/5nu7Y6/wlsPunVWLXOBZ9/VMFM0lGcv0KK2RhOtjyqmEbe\nBJLL64AINkRsnLvafzUg7yHsz3ba77JvTlRXJZD4htTadQPJtTacTBqZVX8KYiA6u7N07k66m27L\nSU6pocJ908jgbYnqNmYSabSD8w/Z5SgY1sspjNbO7qRXDFpeDSjSR15cOdUapnPvLeDuj+PxsQeu\niUDyLYaTGVj7bOWw5tydYmvnitaG03xpu/Lzc0tP7WKYDk0xk6QfgZQSCyepNHOTY/tJrwWV02Ia\nVV+njsIIOyjOJN6f+yCQckXDqVIehdcnrd1Elwgn9wxibxvWwP2wHAuA2mhq2KpWK3lgJZ59W7oD\naqxkkfPQvAmxp7LFLXmX7oAKaatG5zm1SKN0uWMT8XO1e7X6L/Z8+JBfnDcSLiq/nfBG3Q2BVFOV\ncNo4aVQwu+N10xXs0V6ADzm6rrv//df8Qvtzo3cI3Q1qEUitrAqnT3WSXRjVmt1J9y+0CCeG6SCR\nNBIR98IdwgmqEEh7yAyneRo13XRVVvZ2VwknCiOYj0AwjTztwokiSSeaGr5pPoDuafqBqbUEQ344\nkUaIFUar3P/zn/bn4ndUZiDxpt0TgaRFu1XJ5+ouwTC/Z2NeUfFmu7gqaeTZEk45mUQg7YlA0qhp\nOO2wNp19IPszb7OLyx+m26IgnBYziUDaE3NIGrXbMqNdGkWnu0ijy2tRGHncKLIPan8ufgeSRjsj\nkLQLhlPxvqsVt9ebiDRZ8HlGxTSap854+d9/uQ9x/89/eo9oPzjzdyPdDaoQSEdiP04F+2Ws7YUL\nXLqmtY9JIxQP0yWCJ+e286vZS4b//W9757w5FWIO6fDyT3KaXx6+x80rBlEYYbEw2pI6iTtcdXP3\nOVDia0AgnYoXTol5nXbr1PEZhptG1YMnfG+hwmjVze1mzZ+dZXgn745AOq09+8gNhungnVrXupdh\nfJRtaWTv5M///ZeIPP75r9+LeDPvjUC6hB3Cia+TiBVGTWOpShqJE0jW4z/+5f4nb+8dEEiX0yKc\nSCMkJo0ahVPBpNHiHdpMevzHvxKjjrzbG6HL7nLqnuTEMB0Wu+m8nuzg5WvVKow8n/E656maB3Lz\nj1VfG6FCwqggnCiMsOU0I++Iv+qG2x/Uv/w3dcaZJKdCSj9ilYX1YBBICMgJJ9IILU56zWkWzzwh\nKXz50m1NJnlDdpl/JuG0EYGEBfNwYpgOTdemi4XTPBha9JS7mbdlmopwKkAgYQVzGOI9c3E7rE1n\nxSJn/G2zRLz//VfFaSrCKROBBGCFdmmUzp7JNdvHYaOmCSGckggkALmqb683uTy4o3lwbbqWJzlV\n7yZffCyDQ7EQSAByrJ00qjK7k1Om1A2nfQqjyeW/Y4OEkxBIABZlnvQ6uXzbYb2sTNkYTrWWIApf\nHqn/gr+6bDgRSABSvOXpXE17Crbfyefn0n0rcu58cnnt+k8uFk4EEoAwN4r2bKtrNC00/hwpRxIV\nTODyL40Nnj6cCCQAAfNhuuKFFTLt003ghVOitU/5n3nKcCKQAPjyt9ereKbOnlEU+O1xdsoI3q39\n+dCHdAIJwMeWbrot68u1OAV1cnmo2mux6mv6ue1zBpX9+XCHdwIJwGjjaUYFx/Ti8astszuLwVA9\nnHYYjQy+IIc7vBNIAERqL8GQc0xfFQyTy3fct2J7OO1f/5lrHvHYTiABV9d0pVSJHNPdw3S7Zrb5\n09i/BNy//pMD1kYGgQRc2m4rpe6TOolHbxFv489rTnoN3sPk8m2ReehDOoEEXFf1NFo8yM43d/Cu\n0ELrboL5H2Ifcc8YPnoaCYEEXNPGYbqy42zmvqvVp7J27m2bXL5j9XmCgzmBBFxOZmFU9yBbtu/q\nnv0LOXcYvjxU/8leJeAJCiOLQAKuJbEEg6dub9j2fVdbrBS3+KCTy0vrP2lZAp7pGE4gAVcRXCZ1\nh/URWjQ9Lx76V+1bMbm8fW93lXUuTjNM5yKQgEswhdFufQTjQ+yyvZC34MLO9V/wQVfdtuD5nKww\nsggk4PyCk0an3HfV/9Uu6+NVLwEXxwbPetwmkICTy2lhOPq+q+Nv9y0Bv1L/yemG6VwEEnBaZb3d\nX993de0gW+wRW5eA+9d/pz9cE0jAOVU56XXVMX3t+NW39l2tkpetd5GYXP73XycepnMRSMAJtVgQ\naHHf1f03XS27ny3h9JX6Ty5QGxkEEnAqrVdKNTTsu1pxYwjJe55lvd3+5XTTxRFIwHm0Xik13VMg\nO7bVNR0xi3UkpjsVJ5dXCsurHZ8JJOAkKqZR/r6riZN+lNdGi48y/vCl+k8uM0znIpCAw9veTTe5\nvNJOP4fbd1VP/XfZwzKBBBzbYmFUfVipIBiOuO/q/M/cp/678jGZQAIOzKbRbtsf1G0za7fv6uHq\nP7l8GgmBBByXu1jqQbsJ1O67un/9J5ecNPIQSMCx7RNLO6xNp2Hf1e1/Zln9x3HYIJCA82gUTl9Z\nJmdy+UnrPyGNpggk4JxqhZOefVdb12f713/CMN0UgQScX3E4qd13tW79V/cO0w9kceydI5CAa8kM\np+vsu7pz/ScM08URSMB1xcJp/31Xt5QpxU1uX6n/hNoojkACIDINJ8+BdrrLC1ZjCAAAAx9JREFU\nD6ecB61e/3G8TSOQAPj2PMPpW/uuxpZgmNxD1ZkqDraLCCQAKe3CabduAolHzg5PgGG6fAQSgFw2\nnKpsr7dDYTS5/O+/zIPusCSd+2Q4xuYjkACstrFs+ta+q+kNjVqs+8ABdhUCCcAmq8Jp7TBdrX1X\nMzeELbv/2F1xdF2LQAJQTTqcYsHQet/Vshqu+GlQGBUjkAA04YXTt/ZdrTLd9fk5owTkoFqMQALQ\n3Al2yrB3+/k5dO4wR9QtCCQAu9ohnPZZnpy16aojkAB8TfVwan1uU2zgkQNpFQQSABW2h1OjJYg+\nF0a6xjmK1kIgAVCnIJwK0mhjdx/9C9URSABUWwyn9DBdo55y0qgFAgnAYczDyRZGe25/Tho1QiAB\nOKT9W8mFSaPGCCQAh7dPOFEYtUYgATiVRuFEGu2AQAJwWlXCiWG63RBIAC6hLJwojPZEIAG4nMxw\nIo12RiABuLRgODFM9xUEEgCMbDhxYPwKAgkAoMK/ffsJAAAgQiABAJQgkAAAKhBIAAAVCCQAgAoE\nEgBABQIJAKACgQQAUIFAAgCoQCABAFQgkAAAKhBIAAAVCCQAgAoEEgBABQIJAKACgQQAUIFAAgCo\nQCABAFQgkAAAKhBIAAAVCCQAgAoEEgBABQIJAKACgQQAUIFAAgCoQCABAFQgkAAAKhBIAAAVCCQA\ngAoEEgBABQIJAKACgQQAUIFAAgCoQCABAFQgkAAAKhBIAAAVCCQAgAoEEgBABQIJAKACgQQAUIFA\nAgCoQCABAFQgkAAAKhBIAAAVCCQAgAoEEgBABQIJAKACgQQAUIFAAgCoQCABAFQgkAAAKhBIAAAV\nCCQAgAoEEgBABQIJAKACgQQAUIFAAgCoQCABAFQgkAAAKhBIAAAVCCQAgAoEEgBABQIJAKACgQQA\nUIFAAgCoQCABAFQgkAAAKhBIAAAVCCQAgAoEEgBABQIJAKACgQQAUIFAAgCoQCABAFQgkAAAKhBI\nAAAVCCQAgAoEEgBABQIJAKACgQQAUIFAAgCoQCABAFQgkAAAKhBIAAAVCCQAgAoEEgBAhf8H/hU1\n82DDcr4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "% 2-D hat basis\n",
    "clf; \n",
    "[node,elem] = squaremesh([-1,1,-1,1],0.5);\n",
    "phi = zeros(size(node,1),1);\n",
    "phi(12) = 1;\n",
    "subplot(1,2,1);\n",
    "showmesh(node,elem,'facecolor','none'); hold on;\n",
    "showsolution(node,elem,phi,[30,26],'facecolor','g','facealpha',0.5,'edgecolor','k');\n",
    "u = 3-node(:,1).^2 - node(:,2).^2; % x^2 + y^2\n",
    "% u = sin(node(:,1)).*cos(node(:,2));\n",
    "subplot(1,2,2);\n",
    "showmesh(node,elem); hold on;\n",
    "showsolution(node,elem,u,[30,26],'facecolor','g','facealpha',0.5,'edgecolor','k');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Linear Algebraic System\n",
    "\n",
    "For any function $v\\in V_{h}$, there is a unique representation: $v=\\sum _{i=1}^Nv_i\\phi _i$.  We thus can define an isomorphism $V_{h}\\cong \\mathbb R^N$ by\n",
    "\n",
    "$$\n",
    "v=\\sum _{i=1}^Nv_i\\phi _i \\longleftrightarrow \\boldsymbol  v=(v_1, \\cdots, v_N)^{\\intercal},\n",
    "$$\n",
    "\n",
    "and call $\\boldsymbol  v$ the coordinate vector of $v$ relative to the basis $\\{\\phi _i\\}_{i=1}^{N}$. Following the terminology in elasticity, we introduce the *stiffness matrix*\n",
    "\n",
    "$$\n",
    "\\boldsymbol  A=(a_{ij})_{N\\times N}, \\, \\text{ with } \\quad a_{ij}=a(\\phi _j,\\phi _i),\n",
    "$$\n",
    "\n",
    "and the load vector $\\boldsymbol  f=\\{\\langle f, \\phi_k \\rangle\\}_{k=1}^{N}\\in \\mathbb{R}^{N}$. Then the coefficient vector can be obtained by solving the following linear algebraic system\n",
    "\n",
    "$$\n",
    "\\boldsymbol  A\\boldsymbol  u = \\boldsymbol  f.\n",
    "$$\n",
    "\n",
    "It is straightforward to verify $\\boldsymbol  A$ is a symmetric and positive definite (SPD) matrix and thus the solution $\\boldsymbol  u$ exists and unique."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Further reading\n",
    "\n",
    "- Error estimate of the linear finite element method can be found in [Introduction to Finite Element Methods](http://www.math.uci.edu/~chenlong/226/Ch2FEM.pdf)\n",
    "\n",
    "- Assembling the matrix $\\boldsymbol A$ and the vector $\\boldsymbol f$ can be found in [Progamming of Finite Element Methods](http://www.math.uci.edu/~chenlong/226/Ch3FEMCode.pdf)\n",
    "\n",
    "- Fast solvers for solving the linear algebraic equation can be found in [Multigrid Methods](https://www.math.uci.edu/~chenlong/226/MGintroduction.pdf)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Matlab",
   "language": "matlab",
   "name": "matlab"
  },
  "language_info": {
   "codemirror_mode": "octave",
   "file_extension": ".m",
   "help_links": [
    {
     "text": "MetaKernel Magics",
     "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
    }
   ],
   "mimetype": "text/x-matlab",
   "name": "matlab",
   "version": "0.14.3"
  },
  "toc": {
   "colors": {
    "hover_highlight": "#DAA520",
    "navigate_num": "#000000",
    "navigate_text": "#333333",
    "running_highlight": "#FF0000",
    "selected_highlight": "#FFD700",
    "sidebar_border": "#EEEEEE",
    "wrapper_background": "#FFFFFF"
   },
   "moveMenuLeft": true,
   "nav_menu": {
    "height": "102px",
    "width": "252px"
   },
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 4,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": true,
   "widenNotebook": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
